Лабораторная работа 3¶

Сыромятников Дмитрий (КН-301)

In [1]:
N = 15

Задание¶

Решить линейную краевую задачу

$$ \begin{cases} y''=y + 2α + 2αx(1-x), x \in [0, 1]\\ y'(0) = -α\\ y'(1) = e - \frac{1}{e} + α \end{cases} $$

  1. Аналитически
  2. Методом стрельбы
  3. Методом разностной прогонки
  • Для каждого метода строить сетку из N = (5, 10, 50, 100) узлов
  • В решениях методом стрельбы:
    • Сведенную к задаче коши систему уравнений при подброре параметра λ вычислять методами:
      • Эйлера (явный)
      • Эйлера с пересчетом
      • Рунге-Кутта 4-го порядка
    • Корень уравнения $\phi(y'(1|λ) = 0$, где $ϕ$ - функция невязки искать с погрешностью $ϵ=10^{-5}$
    • Очередное значение параметра λ пересчитывать, используя метод Ньютона. Производную невязки приближенно вычислять, используя формулу центральной разностной аппроксимации с приращением $δ=10^{-6}$
  • В решениях методом разностной прогонки проводить аппроксимацию производных:
    • Во внутренних узлах через центральные разности
    • В крайних узлах:
      • Через прямую разность первого порядка
      • Через центральные разности, введя фиктивные узлы

Для каждого N построить графики решений в одной системе координат. Сделать выводы о скорости сходимости методов с ростом числа узлов.

Аналитическое решение¶

$y'' - y = 2α + 2αx(1-x)$

Общее однородное решение¶

$ y'' - y = 0\\ λ^2 - 1 = 0\\ λ = \pm 1\\ ⇒ y_{оо} = C_1e^{x} + C_2e^{-x} $

Частное неоднородное решение¶

$ y_{чн} = px^2 + qx + r\\ y'_{чн} = 2px + q\\ y''_{чн} = 2p\\ 2p - px^2 - qx - r = 2α + 2αx(1-x) = -2α x^2 + 2αx + 2α\\ -px^2-qx+(2p-r) = -2αx^2+2ax + 2α\\ ⇒ \begin{cases} -p = -2α\\ -q = 2α\\ 2p-r = 2α \end{cases} $

$⇒p = 2α, q=-2α, r=2α$ $⇒ y_{чн} = 2αx^2-2αx+2α$

Решение с параметрами¶

$ ⇒ y_{он} = C_1e^{x} + C_2e^{-x} + 2αx^2 - 2αx + 2α $

Поиск параметров¶

$y'_{он} = C_1e^{x} - C_2e^{-x} + 4αx - 2α\\ \begin{cases} y'_{он}(0) = -α ⇒ C_1 - C_2 = α ⇒ C_1 = C_2 + α\\ y'_{он}(1) = e - \frac{1}{e} + \alpha ⇒ C_1e - \frac {C_2}{e} = e - \frac{1}{e} - \alpha ⇒ C_2(e - \frac{1}{e}) = -α(e+1) + e - \frac{1}{e} ⇒ C_2 = -α\frac{e+1} {e - \frac{1}{e}} + 1 \end{cases} $ $ C_1 = 1 - α\frac{e+1} {e - \frac{1}{e}} + α\\ C_2 = 1 - α\frac{e+1} {e - \frac{1}{e}} $

Искомое решение¶

$y = (1 - α(\frac{e+1} {e - \frac{1}{e}} - 1))e^x + (1 - α\frac{e+1} {e - \frac{1}{e}})e^{-x} + 2α(x^2 - x + 1)$

In [2]:
from math import exp, e
In [3]:
α = 2 + 0.1 * N
a, b = 0, 1
ϵ = 1e-5
Ns = (5, 10, 50, 100)
δ = 1e-6
In [4]:
def y_exact(x):
    """
    Точное решение краевой задачи
    y'' = y + 2α + 2α x (1 - x)
    с условиями y'(0) = -α, y'(1) = e - 1/e + α

    :param x: точка, в которой вычисляется решение
    :return: значение y(x)
    """
    tmp = (e + 1) / (e - (1 / e))
    C1 = 1 - α * (tmp - 1)
    C2 = 1 - α * tmp

    return C1 * exp(x) + C2 * exp(-x) + 2 * α * (x*x - x + 1)
In [5]:
x_values = [a + i * (b - a) / 999 for i in range(1000)]
y_exact_values = [y_exact(x) for x in x_values]

Метод стрельбы¶

Сведение задачи к семейству задач коши с параметром λ¶

Сделаем замену $u=y'$. Тогда данную краевую задачу можно переписать в виде задачи коши c параметром λ: $$ \begin{cases} y' = u\\ u' = y + 2α + 2αx(1-x)\\ y(0) = λ\\ u(0) = -α\\ \end{cases} $$

Функция невязки¶

Условие $y'(1) = e - \frac{1}{e} + \alpha$, используем для определения функции невязки, чтобы найти параметр λ:

$ϕ(λ) = y'(1|λ) - y'(1) = y'(1|λ) - (e - \frac{1}{e} + \alpha)$

где $y'(1|λ)$ - приближенное значение производной искомой функции в 1, вычисленное при решении системы с зафиксированным параметром λ

Методы решения задачи Коши системы 2 дифференциальных уравнений¶

Перечисление функций, решающих полученную однопараметризованную систему требуемыми методами:

Явный метод Эйлера¶

In [6]:
def euler_explicit(f, a, b, y0, u0, n):
    """
    Численно решает задачу Коши для системы из двух обыкновенных дифференциальных уравнений:
        y' = u,
        u' = f(x, y)

    Используется явный метод Эйлера (метод прямого шага по x) на равномерной сетке из n узлов.

    Параметры:
        f  : функция f(x, y), задающая правую часть второго уравнения u' = f(x, y)
        a  : float, левая граница интервала интегрирования
        b  : float, правая граница интервала интегрирования
        y0 : float, начальное значение функции y в точке x = a
        u0 : float, начальное значение функции u в точке x = a
        n  : int, количество узлов сетки (в том числе начальный и конечный)

    Возвращает:
        X : list of float, координаты узлов сетки на интервале [a, b]
        U : list of float, приближённые значения функции u(x) = y'(x) в узлах X
        Y : list of float, приближённые значения функции y(x) в узлах X
    """
    h = (b - a) / (n - 1)
    X = [a + i * h for i in range(n)]
    U = [u0]
    Y = [y0]
    for i in range(n - 1):
        U.append(U[i] + h * f(X[i], Y[i]))
        Y.append(Y[i] + h * U[i])
    return X, U, Y

Метод Эйлера с пересчетом¶

In [7]:
def euler_recalc(f, a, b, y0, u0, n):
    """
    Численно решает задачу Коши для системы из двух обыкновенных дифференциальных уравнений:
        y' = u,
        u' = f(x, y)

    Используется метод Эйлера c пересчетом на равномерной сетке из n узлов.

    Параметры:
        f  : функция f(x, y), задающая правую часть второго уравнения u' = f(x, y)
        a  : float, левая граница интервала интегрирования
        b  : float, правая граница интервала интегрирования
        y0 : float, начальное значение функции y в точке x = a
        u0 : float, начальное значение функции u в точке x = a
        n  : int, количество узлов сетки (в том числе начальный и конечный)

    Возвращает:
        X : list of float, координаты узлов сетки на интервале [a, b]
        U : list of float, приближённые значения функции u(x) = y'(x) в узлах X
        Y : list of float, приближённые значения функции y(x) в узлах X
    """
    h = (b - a) / (n - 1)
    X = [a + i*h for i in range(n)]
    U = [u0]
    Y = [y0]
    for i in range(n - 1):
        xi, yi, ui = X[i], Y[i], U[i]
        fi = f(xi, yi)

        us = ui + h * fi
        ys = yi + h * ui

        fs = f(X[i+1], ys)
        u_next = ui + (h / 2) * (fi + fs)
        y_next = yi + (h / 2) * (ui + us)
        U.append(u_next)
        Y.append(y_next)
    return X, U, Y

Метод Рунге-Кутта 4-го порядка¶

In [8]:
def runge_kutta(f, a, b, y0, u0, n):
    """
    Численно решает задачу Коши для системы из двух обыкновенных дифференциальных уравнений:
        y' = u,
        u' = f(x, y)

    Используется метод Рунге-Кутта 4-го порядка на равномерной сетке из n узлов.

    Параметры:
        f  : функция f(x, y), задающая правую часть второго уравнения u' = f(x, y)
        a  : float, левая граница интервала интегрирования
        b  : float, правая граница интервала интегрирования
        y0 : float, начальное значение функции y в точке x = a
        u0 : float, начальное значение функции u в точке x = a
        n  : int, количество узлов сетки (в том числе начальный и конечный)

    Возвращает:
        X : list of float, координаты узлов сетки на интервале [a, b]
        U : list of float, приближённые значения функции u(x) = y'(x) в узлах X
        Y : list of float, приближённые значения функции y(x) в узлах X
    """
    h = (b - a) / (n - 1)

    X = [a + i * h for i in range(n)]
    Y = [y0]
    U = [u0]

    for i in range(n - 1):
        x = X[i]
        y = Y[i]
        u = U[i]

        k1y = u
        k1u = f(x, y)

        k2y = u + 0.5 * h * k1u
        k2u = f(x + 0.5 * h, y + 0.5 * h * k1y)

        k3y = u + 0.5 * h * k2u
        k3u = f(x + 0.5 * h, y + 0.5 * h * k2y)

        k4y = u + h * k3u
        k4u = f(x + h, y + h * k3y)

        y_next = y + (h / 6) * (k1y + 2*k2y + 2*k3y + k4y)
        u_next = u + (h / 6) * (k1u + 2*k2u + 2*k3u + k4u)

        Y.append(y_next)
        U.append(u_next)

    return X, U, Y

Подбор параметра λ¶

Начальным $λ$ возьмем $0$. Будем перерешивать систему и уточнять λ, пересчитывая его с использованием метода Ньютона, пока функция невязки $ϕ$ от $λ$ не станет близка к $0$ (корню) с требуемой точностью $ϵ$:

In [9]:
def f(x, y):
  """
  функция f(x, y), задающая правую часть второго уравнения
  """
  return y + 2*α + 2*α*x*(1-x)
In [10]:
def ϕ(λ, method):
    """
    Функция невязки и метод, с помощью которого она построена
    """
    X, U, Y = method(f, a, b, λ, -α, n)
    return U[-1] - (e - 1/e + α)

def dϕ(λ, method):
    """
    Производная функции невязки и метод, с помощью которого она будет найдена
    """
    return (ϕ(λ + δ, method) - ϕ(λ - δ, method)) / (2 * δ)

Явный метод Эйлера¶

In [11]:
X_euler_s = []
Y_euler_s = []
for n in Ns:
  λ = 0
  while True:
      ϕλ = ϕ(λ, euler_explicit)
      if abs(ϕλ) < ϵ:
          break
      dϕλ = dϕ(λ, euler_explicit)
      if dϕλ == 0:
          raise ZeroDivisionError()

      λ = λ - ϕλ / dϕλ

  X_euler, _, Y_euler = euler_explicit(f, a, b, λ, -α, n)
  X_euler_s.append(X_euler)
  Y_euler_s.append(Y_euler)

Метод Эйлера с пересчетом¶

In [12]:
X_euler_recalc_s = []
Y_euler_recalc_s = []
for n in Ns:
  λ = 0
  while True:
      ϕλ = ϕ(λ, euler_recalc)
      if abs(ϕλ) < ϵ:
          break
      dϕλ = dϕ(λ, euler_recalc)
      if dϕλ == 0:
          raise ZeroDivisionError()

      λ = λ - ϕλ / dϕλ

  X_euler_recalc, _, Y_euler_recalc = euler_recalc(f, a, b, λ, -α, n)
  X_euler_recalc_s.append(X_euler_recalc)
  Y_euler_recalc_s.append(Y_euler_recalc)

Метод Рунге-Кутта 4-го порядка¶

In [13]:
X_runge_kutta_s = []
Y_runge_kutta_s = []
for n in Ns:
  λ = 0
  while True:
      ϕλ = ϕ(λ, runge_kutta)
      if abs(ϕλ) < ϵ:
          break
      dϕλ = dϕ(λ, runge_kutta)
      if dϕλ == 0:
          raise ZeroDivisionError()

      λ = λ - ϕλ / dϕλ

  X_runge_kutta, _, Y_runge_kutta = runge_kutta(f, a, b, λ, -α, n)
  X_runge_kutta_s.append(X_runge_kutta)
  Y_runge_kutta_s.append(Y_runge_kutta)

Метод разностной прогонки¶

Дискретезация области решения¶

Разделим отрезок $[0, 1]$ на $N$ равных частей с шагом $h=\frac {1}{N}$

Получим узлы: $x_0 = 0, x_1 = h, x_2 = 2h, ..., x_N = 1$

Введем обозначение: $y_i ≈ y(x_i)$ - приближенное значение функции в узле

Аппроксимация производных во внутренних узлах рассматриваемого отрезка¶

Для всех внутренних узлов: $x_1, x_2, ... x_{Ns - 1}$

Апроксимируем $y''(x_i)$ через формулы центральных разностей: $$ y''(x_i) ≈ \frac{y_{i-1} - 2y_i + y_{i+1}}{h^2}, i = 1...(N-1) $$

Подставим в исходное уравнение: $$ \frac{y_{i-1} - 2y_i + y_{i+1}}{h^2} = y_i + 2α + 2αx_i(1-x_i) $$

Домножим на $h^2$, перенесем влево $h^2y_i$, приведем подобные. Получим разностное уравнение: $$ y_{i-1} -(2+h^2)y_i + y_{i+1} = 2αh^2(1 + x_i(1-x_i)), i = 1...(N-1) $$

Получили линейную относительно $y_i$ систему. Подробнее:

$$ \begin{cases} y_0 -(2+h^2)y_1 + y_2 = 2αh^2(1 + x_1(1-x_1))\\ y_1 -(2+h^2)y_2 + y_3 = 2αh^2(1 + x_2(1-x_2))\\ ...\\ y_{i-1} -(2+h^2)y_i + y_{i+1} = 2αh^2(1 + x_i(1-x_i))\\ ...\\ y_{N-2} -(2+h^2)y_{N-1} + y_N = 2αh^2(1 + x_{N-1}(1-x_{N-1})) \end{cases} $$

Аппроксимация краев односторонними разностями¶

В $x_0$: $$ y'(0) ≈ \frac {y_1 - y_0} {h} = -α ⇒ y_1 - y_0 = -αh $$

В $x_N$: $$ y'(1) ≈ \frac {y_N - y_{N-1}} {h} = e - \frac{1}{e} + α ⇒ y_N - y_{N-1} = h(e - \frac{1}{e} + α) $$

Формирование системы¶

Объединяя систему, построенную для внутренних узлов с уравнениями, выведенными из краевых условий c использованием односторонних разностей, получаем систему из $N+1$ уравнения с $N+1$ неизвестными $y_i , i = 0,...N$:

$$ \begin{cases} y_1 - y_0 = -αh\\ y_0 -(2+h^2)y_1 + y_2 = 2αh^2(1 + x_1(1-x_1))\\ y_1 -(2+h^2)y_2 + y_3 = 2αh^2(1 + x_2(1-x_2))\\ ...\\ y_{i-1} -(2+h^2)y_i + y_{i+1} = 2αh^2(1 + x_i(1-x_i))\\ ...\\ y_{N-2} -(2+h^2)y_{N-1} + y_N = 2αh^2(1 + x_{N-1}(1-x_{N-1}))\\ y_N - y_{N-1} = h(e - \frac{1}{e} + α) \end{cases} $$

Обозначим $ξ = -(2+h^2), χ =2αh^2$

Тогда в матричном представлении данная система выглядит следующим образом: $$ \begin{bmatrix} -1 & 1 & 0 & 0 & \cdots & 0 & 0 \\ 1 & ξ & 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & ξ & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & ξ & 1 & 0 \\ 0 & 0 & 0 & \cdots & 1 & ξ & 1 \\ 0 & 0 & 0 & \cdots & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-2} \\ y_{N-1} \\ y_N \end{bmatrix} = \begin{bmatrix} -\alpha h \\ χ\bigl(1 + x_1(1-x_1)\bigr) \\ χ\bigl(1 + x_2(1-x_2)\bigr) \\ \vdots \\ χ\bigl(1 + x_{N-2}(1-x_{N-2})\bigr) \\ χ\bigl(1 + x_{N-1}(1-x_{N-1})\bigr) \\ h\bigl(e - \tfrac{1}{e} + \alpha\bigr) \end{bmatrix} $$

Решение системы¶

Матрица системы трехдиагональная. Поэтому можем решить эту систему методом прямой и обратной прогонки для всех заданных N

In [14]:
def thomas(As, Bs, Cs, Ds):
    """
    Решает систему линейных уравнений Ax = D, где A — трехдиагональная матрица,
    заданная в виде трех списков чисел:
      - As — нижняя диагональ (считается, что A[0] = 0),
      - Bs — главная диагональ,
      - Cs — верхняя диагональ (считается, что C[n-1] = 0),
    методом Томаса (прямой и обратной прогонки).

    Параметры:
        As: нижняя диагональ (размер n)
        Bs: главная диагональ (размер n)
        Cs: верхняя диагональ (размер n)
        Ds: правая часть системы (размер n)

    Возвращает:
        решение системы линейных уравнений - вектор Y размера n
    """
    Ps = [0] * n
    Qs = [0] * n
    Y  = [0] * n
    # Прямая прогонка
    Ps[0] = -Cs[0] / Bs[0]
    Qs[0] =  Ds[0] / Bs[0]

    for i in range(1, n):
        denom = As[i] * Ps[i-1] + Bs[i]
        if i < n-1:
            Ps[i] = -Cs[i] / denom
        Qs[i] = (Ds[i] - As[i] * Qs[i-1]) / denom

    # Обратная прогонка
    Y[n-1] = Qs[n-1]
    for i in range(n-2, -1, -1):
        Y[i] = Ps[i] * Y[i+1] + Qs[i]
    return Y
In [15]:
race = []

for n in Ns:
    # 1. Сетка
    h = (b - a) / (n - 1)
    X = [a + i * h for i in range(n)]

    ξ = -(2 + h * h)
    χ = 2 * α * h * h

    As = [0] + [1] * (n - 2) + [-1]
    Bs = [-1] + [ξ] * (n - 2) + [1]
    Cs = [1] * (n - 1) + [0]

    Ds = [ -α * h ]
    for xi in X[1:-1]:
        Ds.append( χ * (1 + xi * (1 - xi)))
    Ds.append( h * (e - 1/e + α) )

    race.append((X, thomas(As, Bs, Cs, Ds)))

Аппроксимация через введение фиктивных узлов¶

Аппроксимируем значения производной искомой функции на концах рассматриваемого отрезка через формулы центральных разностей, введя фиктивные узлы, лежащае за границами рассматриваемого отрезка $[a, b]$: $x_{-1} = a - h, x_{N+1} = b + h$ (формально)

В $x_0$:

$$ y'(x_0) ≈ \frac {y_{-1} - 2y_0 + y_1}{h^2} $$

В $x_N$: $$ y'(x_N) ≈ \frac {y_{N-1} - 2y_N + y_{N+1}}{h^2} $$

Приблеженное значение функции в фиктивных узлах выразим из краевых условий и формулы приближения производной первого порядка (центральной):

$y_{-1}$: $$ y'(x_0) ≈ \frac {y_1 - y_{-1}} {2h} = -α ⇒ y_{-1} = y_1 + 2αh $$

$y_{N+1}$: $$ y'(x_N) ≈ \frac {y_{N+1} - y_{N-1}} {2h} = e - \frac{1}{e} + α ⇒ y_{N+1} = y_{N-1} + 2h(e - \frac{1}{e} + α) $$

Подставим выраженные значения в апроксимированные крайние производные, упростим. Приравнивая полученное выражение к данному в условии уравнению $y''$, получим уравнения:

$x_0$: $$ \frac {y_1 + 2αh - 2y_0 + y_1}{h^2} = y_0 + 2α + 2αx_0(1-x_0) ⇒ -(h^2 + 2)y_0 + 2y_1 = 2α(h^2(1 + x_0(1-x_0)) - h) $$

$x_N$: $$ \frac {y_{N-1} - 2y_N + y_{N-1} + 2h(e - \frac{1}{e} + α)}{h^2} = y_N + 2α + 2αx_N(1-x_N) ⇒ 2y_{N-1} -(h^2 + 2)y_N = 2αh^2(1 + x_N(1-x_N)) - 2h(e - \frac{1}{e} + α) $$

Формирование системы¶

Объединяя систему, построенную для внутренних узлов с уравнениями, выведенными выше, получаем систему из $N+1$ уравнения с $N+1$ неизвестными $y_i , i = 0,...N$:

$$ \begin{cases} -(2 + h^2)y_0 + 2y_1 = 2α(h^2(1 + x_0(1-x_0)) - h) y_0 -(2+h^2)y_1 + y_2 = 2αh^2(1 + x_1(1-x_1))\\ y_1 -(2+h^2)y_2 + y_3 = 2αh^2(1 + x_2(1-x_2))\\ ...\\ y_{i-1} -(2+h^2)y_i + y_{i+1} = 2αh^2(1 + x_i(1-x_i))\\ ...\\ y_{N-2} -(2+h^2)y_{N-1} + y_N = 2αh^2(1 + x_{N-1}(1-x_{N-1}))\\ 2y_{N-1} -(2 + h^2)y_N = 2αh^2(1 + x_N(1-x_N)) - 2h(e - \frac{1}{e} + α) \end{cases} $$

Обозначим $ξ = -(2+h^2), χ =2αh^2$

Тогда в матричном представлении данная система выглядит следующим образом:

$$ \begin{bmatrix} \xi & 2 & 0 & 0 & \cdots & 0 & 0\\ 1 & \xi & 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \xi & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & \xi & 1 & 0 \\ 0 & 0 & 0 & \cdots & 1 & \xi & 2 \\ 0 & 0 & 0 & \cdots & 0 & 2 & \xi \end{bmatrix} \begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{N-1}\\ y_N \end{bmatrix} = \begin{bmatrix} 2\alpha\bigl(h^2(1 + x_0(1-x_0)) - h\bigr)\\[6pt] \chi\bigl(1 + x_1(1-x_1)\bigr)\\[6pt] \chi\bigl(1 + x_2(1-x_2)\bigr)\\[3pt] \vdots\\[3pt] \chi\bigl(1 + x_{N-1}(1-x_{N-1})\bigr)\\[6pt] \chi\bigl(1 + x_N(1-x_N)\bigr) \;-\; 2h\Bigl(e - \tfrac1e + \alpha\Bigr) \end{bmatrix} $$

Решение системы¶

Матрица системы трехдиагональная. Поэтому можем решить эту систему методом прямой и обратной прогонки для всех заданных N

In [16]:
race_ghost = []

for n in Ns:
    h = (b - a) / (n - 1)
    X = [a + i*h for i in range(n)]
    ξ = -(2 + h*h)
    χ =  2*α*h*h

    As = [0] + [1]*(n-2) + [2]
    Bs = [ξ]*n
    Cs = [2] + [1]*(n-2) + [0]

    Ds = [2*α*(h*h*(1 + X[0]*(1 - X[0])) - h)]
    for xi in X[1:-1]:
        Ds.append(χ*(1 + xi*(1 - xi)))
    Ds.append(χ*(1 + X[-1]*(1 - X[-1])) - 2*h*(e - 1/e + α))

    race_ghost.append((X, thomas(As, Bs, Cs, Ds)))

Графики¶

In [17]:
colors = {
    'euler': 'blue',
    'euler_recalc': 'green',
    'runge_kutta': 'red',
    'race': 'orange',
    'race_ghost': 'yellow'
}
In [18]:
import plotly.graph_objects as go

fig = go.Figure()

fig.add_trace(go.Scatter(
    x=x_values, y=y_exact_values,
    name='Точное решение',
    line=dict(color='black', dash='dash', width=2),
    visible=True
))

for i, n in enumerate(Ns):
    visible_flag = (i == 0)

    fig.add_trace(go.Scatter(
        x=X_euler_s[i], y=Y_euler_s[i],
        name='Эйлер явный',
        line=dict(color=colors['euler']),
        mode='lines+markers',
        visible=visible_flag
    ))

    fig.add_trace(go.Scatter(
        x=X_euler_recalc_s[i], y=Y_euler_recalc_s[i],
        name='Эйлер с пересчетом',
        line=dict(color=colors['euler_recalc']),
        mode='lines+markers',
        visible=visible_flag
    ))

    fig.add_trace(go.Scatter(
        x=X_runge_kutta_s[i], y=Y_runge_kutta_s[i],
        name='Рунге-Кутта',
        line=dict(color=colors['runge_kutta']),
        mode='lines+markers',
        visible=visible_flag
    ))

    fig.add_trace(go.Scatter(
        x=race[i][0], y=race[i][1],
        name='Разностной прогонки',
        line=dict(color=colors['race']),
        mode='lines+markers',
        visible=visible_flag
    ))

    fig.add_trace(go.Scatter(
        x=race_ghost[i][0], y=race_ghost[i][1],
        name='Разностной прогонки с фиктивными узлами',
        line=dict(color=colors['race_ghost']),
        mode='lines+markers',
        visible=visible_flag
    ))

steps = []
for i, n in enumerate(Ns):
    visibility = [True]
    for j in range(len(Ns)):
        visibility += [j == i] * 5
    steps.append(dict(
        method="update",
        args=[{"visible": visibility},
              {"title": f"Сравнение методов при n = {n}"}],
        label=str(n)
    ))

sliders = [dict(
    active=0,
    currentvalue={"prefix": "Количество узлов: "},
    pad={"t": 50},
    steps=steps
)]

fig.update_layout(
    title='Сравнение методов решения краевой задачи',
    xaxis_title='x',
    yaxis_title='y',
    sliders=sliders,
    showlegend=True,
    legend=dict(orientation="h", yanchor="bottom", y=1.02)
)

fig.show()

Выводы¶

Методы успешно решают краевую задачу, при этом их точность закономерно возрастает с ростом числа узлов (уменьшением шага). Метод стрельбы с использованием метода Рунге-Кутта 4-го порядка имеет среди рассмотренных методов наибольшую скорость сходимости. На втором месте по скорости сходимости идет метод разностной прогонки с аппроксимацией производных на краях через введение фиктивных узлов. Хуже всех в плане скорости сходимости показал себя метод разностной прогонки с аппроксимацией производных односторонними разностями. Вид графиков решения полностью согласуются с теоретическими выкладками о методах, с помощью которых приближенные решения были найдены.